﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.MinPack
{
    /// <summary>
    /// Computes a forward-difference approximation to the m by n jacobian matrix associated with a 
    /// specified problem of m functions in n variables.
    /// </summary>
    [Serializable]
    public static class Fdjac2Class
    {
        /// <summary>
        /// Computes a forward-difference approximation to the m by n jacobian matrix associated with a 
        /// specified problem of m functions in n variables.
        /// </summary>
        /// <param name="f">The name of the user-supplied subroutine which calculates the functions.</param>
        /// <param name="m">Input variable set to the number of functions.</param>
        /// <param name="n">Input variable set to the number of variables. n must not exceed m.</param>
        /// <param name="x">An input array of length n.</param>
        /// <param name="fvec">An input array of length m which must contain the functions evaluated at x.</param>
        /// <param name="fjac">An output m by n array which contains the approximation to the jacobian matrix 
        /// evaluated at x.</param>
        /// <param name="iflag">An integer variable which can be used to terminate the execution of fdjac2. 
        /// See description of fcn.</param>
        /// <param name="epsfcn">An input variable used in determining a suitable step length for the 
        /// forward-difference approximation.</param>
        /// <param name="wa">A work array of length m.</param>
        public static void Fdjac2(MinPackFunction f, int m,
                                  int n, double[] x, double[] fvec, double[,] fjac,
                                  ref int iflag, double epsfcn, double[] wa)
        {
            double[] dpmpar = {
                                  2.220446049250313e-16,
                                  2.225073858507201e-308,
                                  1.797693134862316e-308
                              };
            int i, j;
            double h, temp;
            double epsmch = (epsfcn > dpmpar[0]) ? epsfcn : dpmpar[0];
            double eps = Math.Sqrt(epsmch);

            for (j = 0; j < n; j++)
            {
                temp = x[j];
                if (temp == 0.0)
                {
                    h = eps;
                }
                else
                {
                    h = eps * Math.Abs(temp);
                }
                x[j] = temp + h;
                f(m, n, x, wa, ref iflag);
                if (iflag < 0)
                {
                    break;
                }
                x[j] = temp;
                for (i = 0; i < m; i++)
                {
                    fjac[j, i] = (wa[i] - fvec[i]) / h;
                }
            }
        }
    }
}